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Abstract. This paper focusses on the barely understood gap between the weakly nonlinear regime of structure 
formation and the onset of the virialized regime. While the former is accessed through perturbative calcula- 
tions and the latter through virialization conditions incorporating dynamical stresses that arise in collisionless 
self-gravitating systems due to velocity dispersion forces, the addressed regime can only be understood through 
non-perturbative models. We here present an exact Lagrangian integral that provides a tool to access this regime. 
We derive a transport equation for the peculiar-gravitational field strength and integrate it along comoving tra- 
jectories of fluid elements. The so-obtained integral provides an exact expression that solves the longitudinal 
gravitational field equation in general. We argue that this integral provides a powerful approximation beyond 
the Lagrangian perturbative regime, and discuss its relation to known approximations, among them Lagrangian 
perturbation solutions including the Zel'dovich approximation and approximations for adhesive gravitational clus- 
tering, including the adhesion approximation. Furthermore, we propose an iteration scheme for a systematic an- 
alytical and numerical construction of trajectory fields. The integral may also be employed to improve inverse 
reconstruction techniques. 
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1. Introduction 



Both the analytical understanding and the numerical 
simulation of the formation of cosmic structure in the 
Universe have substantially advanced in recent years. 
Early results in the development of both approaches and 
their comparison has improved our understanding of the 
building blocks of large-scale structure. Both focussed on 
modelling dark matter in terms of a Newtonian dust con- 
tinuum, but recent developments are heading in some- 
what different directions. Numerical attempts are directed 
towards improving spatial resolution, understanding the 
force distribution on a lattice, and incorporating new 
physics for the modelling of structure formation on galaxy 
halo scales (e.g., Bagla & Padmanabhan 1997; Gabrielli 
et al. 2006; Bertschinger 1998; Shirokov & Bertschinger 
2006) and for the modelling of hydrodynamical effects 
(Steinmetz 2003). Analytical improvements have still fol- 
lowed the route of understanding structure formation in 
a self-gravitating dust continuum, but also attempt to 
access small-scale structure by incorporating additional 
forces in the Euler equation that arise from kinetic the- 
ory (e.g., Ma & Bertschinger 2004; Buchert & Dominguez 
2005). In this work we discuss a regime that is difficult to 
understand and to model using both numerical and ana- 
lytical strategies. The analytical insight will also shed light 



on and provide tools for numerical and semi-analytical 
techniques. 

In Sect. |2 we recall the basic equations, as well as cos- 
mological models that have been obtained previously and 
that help to locate the non-perturbative regime from the 
analytical point of view. In Sect. [3| we then formulate the 
key equation of the present work and integrate this equa- 
tion exactly along the flow lines of continuum elements. 
The implications of this result are investigated in Sect.QJ 
and in Sect. [S] we give a discussion summary. 

In this work we do not provide a comprehensive refer- 
ence list, because the reader may find systematic deriva- 
tions of the equations of Sect.[2]as well as a substantial list 
of references in a recent paper on the current status of ana- 
lytical models (Buchert & Domi'nguez 2005) . We first work 
in Eulerian coordinates that are comoving with a given ref- 
erence Hubble flow, i.e. a homogeneous-isotropic solution 
of the Euler-Newton system of equations, as is widely used 
in cosmology. We denote them by q = x/a(t), where x are 
the Eulerian coordinates, a(t) a solution of Friedmann's 
differential equation, and H :— a/a is Hubble's function. 
Later, we move to the Lagrangian picture of fluid motion 
and represent the comoving trajectory field of continuum 
elements by q = F(X, t), where X are the Lagrangian co- 
ordinates, which index fluid elements. The velocity field is 
split into a Hubble-velocity and a peculiar-velocity, v = 



2 
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Vh + u, the acceleration field into a Hubble-acceleration 
and a peculiar-acceleration g = g# + w, the density 
field into a background density and a density contrast 
g = qh{^ + $)■ The total (Lagrangian) time-derivative is 
denoted by A = d t | x + v • V x = d t | q + u/a- V q , and some- 
times by an overdot. The nabla-operator with respect to 
Lagrangian coordinates is denoted by Vo and sometimes 
by V x . 

2. The non-perturbative regime 

In this section we recall the basic equations in a form 
that is helpful for our purpose and then review the 
known cosmological model equations within the presented 
framework and discuss their limits with regard to non- 
perturbative effects in the non-linear regime. 

2.1. Equations for the non-perturbative regime 

For the basic system of equations, derived from coarse- 
graining the Newtonian kinetic equations for a system of 
N self-gravitating particles, we refer the reader to our re- 
cent paper (Buchert & Domfnguez 2005). We work below 
with hydrodynamical equations that can also be found 
in many cosmology textbooks and review papers (e.g., 
Peebles 1980, Binney & Tremaine 1987, Sahni & Coles 
1995, Ehlers & Buchert 1997, Bernardeau et al. 2002). 

We start with recalling equations that describe the evo- 
lution of the density g, the peculiar-velocity u, and the 
peculiar-acceleration or the peculiar-gravitational field 
strength w, obtained by coarse-graining the kinetic equa- 
tions of N self-gravitating particles in real space and ve- 
locity space and then forming velocity moments to obtain 
a hydrodynamical description in Eulerian space. The ze- 
roth velocity moment is the continuity equation, 



d 1 
— g+3Hg+ -V q u 
dt a 







(la) 



The first velocity moment provides the evolution equation 
for the peculiar-velocity: 



4-u + Hu = w + - f F — -V a • n 



dt 



(lb) 



In the above equation, the peculiar-gravitational field 
strength w is constrained by the Newtonian field equa- 
tions: 



V q X W = 



V q • w = —AnGgnaS 



(lc) 



The force F represents deviations from mean field gravity 
(modelled by Eqs. ifTcl) - ). and the symmetric tensor field 
II represents forces due to velocity dispersion. For F = 0, 
Eq. (|lbl) is known as the Euler-Jeans equation in stel- 
lar system theory. The above set of equations truncates 
the velocity moment hierarchy, so we need models for the 
fields F and n (e.g., as functionals of the other fields) 
to close the system of equations. The simplest truncation 
is to neglect deviations from mean field gravity and ve- 
locity dispersion altogether and study the evolution of a 



dust continuum. More general models for adhesive gravi- 
tational clustering including those forces are investigated 
in (Buchert & Domfnguez 2005). 

In this work we focus on the (in the above framework) 
general evolution equation for the peculiar-gravitational 
field strength (Buchert 1989) 



dt 



w + 2iJw - AitGqhu = TZ ; 



-R:=-[(u-Vq)w-u(Vq-w) + V q xT] , (2) 

which can be obtained from Eq. Qlajl by inserting the field 
equations l|Tc|) and formally integrating the divergence. 
The resulting vector field of integration T can be deter- 
mined as follows. Acting with V q x on Eq. © and sub- 
jecting T to the Coulomb gauge condition, V q • T = 0, one 
finds: 

A q T = AnGaVq x (gu) = AnGa [pV q x u + V q £ x u] 

= (AnGagn — V q ■ w) V q x u + u x A q w . (3) 

That is, T is the vector potential (up to an unimportant 
factor) of the peculiar- current density} := gu: 



AirGaj = -V q x T - V q ~<9 t | q [a 



(4) 



with <j) denoting the scalar peculiar-gravitational poten- 
tial. (Note that Eq. is equivalent to Eq. by employ- 
ing the definitions = <9t| q w + u/a ■ V q w and w =: 
— ^Vq<A.) If we require T to vanish, then Equation J3J) 
shows that, for irrotational flows, the mean flow follows 
the gradient of the density field 3 ). 

We thus have obtained a new set of equations that is 
equivalent to the set (jlallcj) . and we can write it in a form 
that reminds us of Maxwell's equations: 



V q X W 



Vq X T = a t | q W-47TGj 



Vq W 



-Vq-T 



—AirGgnS ; 
. (5) 



For the purpose of constructing models for cosmic 
structure formation, one usually employs the equations 
(|lallbllc|) and approximates the fields F and II. A com- 
monly used approximation is to set F equal to zero and 
"masking" Tl by its isotropic contribution, =: pSij, 
together with a dynamical equation of state p = (3(g) 
(Buchert & Domfnguez 1998). For such a closure approx- 
imation, we obtain the starting equations for a large set 
of approximations that have been discussed in the cosmo- 
logical literature. 

To recall these approximations, let us now move to the 
Lagrangian picture of fluid mechanics, where the Eulerian 

3 If such conditions are imposed on the problem, T becomes 
a harmonic vector function that can be set to zero for peri- 
odic boundary conditions, since harmonic functions are then 
spatially constant and can be set to zero due to the invari- 
ance of the basic equations with respect to spatially constant 
translations. 
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positions are viewed as a field variable depending on the 
initial position vectors X (the Lagrangian coordinates) 
and the time t, q = F(X, t). The trajectory field F rep- 
resents the integral curves of the scaled peculiar-velocity 
field, u = aF. We derive equations for the comoving dis- 
placement field P := F X. First, we formally integrate 
Eq. |PJ with respect to the time and obtain (Buchert & 
Dommguez 2005): 



W 1 t l 

-^ + 4irGQ H aP(X,t) + ^ / dt' a 2 (t')K(X, t' 



(6a) 



with w(X,to)=:W(X) , P(X,to):=0 . 

(6b) 

The integration constant has been chosen such that 
the Lagrangian coordinates coincide with the comoving 
Eulerian coordinates at initial time 4 . Notice that Eq. (j6a|) 
is general and does not depend on the specific models for 
the forces in Eq. I)lb[l. Inserting this general expression 
for w into Eq. Ijlbfl and using u = aP yields an evolution 
equation for the displacement field P, which we simplify - 
for the sake of a more transparent discussion - to the case 
F = and an isotropic velocity dispersion tensor (repre- 
senting the dynamical pressure by an equation of state of 
the form p = {3(g)); this restriction implies: 



— v q -n 



(7) 



Above, we have introduced a density dependent Jeans' 
length as the product of the speed of sound c s and the 
local free-fall time tp: 



t t 2 dp i 

Lj := V ~^Gq ^ s F ' s = T e = 13 {g) ' tF «■ 

. (8) 

For this case we obtain the following evolution equation: 
P + 2HP - AttGqhP = 

(9) 



^^Aqw + \ [ dt' a 2 (t')Tl(X, t') 

a a Jtn 



2.2. A list of known cosmological model equations 

From Eq. JJjJ we are able to infer a list of approximation 
assumptions that have been used for the construction of 
cosmological evolution models: 



4 Alternatively, we could define initial displacements through 
47rG{JH-(to)P(X, to) = w(X, to)- This means either that we 
have assumed quasi-homogeneity, p(X, to) ~ Quito) or, else, 
we have adopted a particular choice of Lagrangian coordinates. 
With this assumption we get rid of physically unimportant con- 
stant terms, which is sometimes useful, but will not be used 
here to avoid confusion with the exact property of the expres- 
sions derived in this paper. (For details on the possibility and 
advantage of this choice, see Adler & Buchert 1999, App. A.) 



1 Setting Lj — formally we obtain the dust model with- 
out any restriction. 

2 Neglecting the residual vector field 1Z = 0, we obtain 
for Lj = the equation for the longitudinal part of 
Lagrangian first-order perturbations (the transversal 
part is hidden in 7L = 0). The peculiar-gravitational 
field strength is then, in view of Eq. I|6a|) , approximated 
through the relation: 



W 



AuGqhoP 



(10) 



3 Also neglecting the residual vector field 1Z = 0, we 
obtain for Lj ^ together with the valid relation (|10|) : 



3 A 1 W = 



Qh 13' (q) 



a 2 g 



A q P 



(11) 



The coefficient in front of the (Eulerian) Laplacian of P 
is density-independent for the special choice p = kq 2 , 
k = const., which corresponds to the assumption of 
an isolated and 'virialized' fluid element, cf. Buchert 
& Dommguez 2005, Appendix C). For this example we 
obtain an evolution equation of the form: 

P + 2HP - ^Gq h P = ^|^A q P . (12) 

Note that the above model equation is nonlin- 
ear in Eulerian space (due to the convective non- 
linearities hidden in the overdot) and in Lagrangian 
space (the Eulerian Laplacian, if transformed to 
Lagrangian coordinates is non-linear). It is of a hybrid 
Lagrangian/Eulerian type and has been suggested in 
Buchert & Dommguez (2005). 

4 Linearization of the convective non-linearities in 
Eq. i|12|) yields the Eulerian linear approximation (e.g., 
Peebles 1980). For this, note that linearization of the 
general integral for the density (see Eq. (|20|) below) 
gives S = — Vo ■ P. 

5 By linearizing the Eulerian Laplacian in Eq. (|12(l , now 
in the Lagrangian frame, i.e., retaining only the zero- 
order Lagrangian term q = X + P X, we recover 
the Lagrangian linear approximation for a medium 
supported by a dynamical pressure (Adler & Buchert 
1999). 

6 Restricting the latter model further by assuming that 
u oc w, with the function of proportionality taken from 
the Eulerian linear approximation for a dust contin- 
uum, we obtain the standard adhesion approximation 
(Gurbatov et al. 1989). 

7 The adhesion approximation reduces, for Lj = 0, to 
the ZeVdovich approximation (Zel'dovich 1970, 1973). 

2.3. Discussion of models 

While models 2 and 4-7 are perturbative, model 3 already 
describes the non-perturbative regime concerning the dy- 
namical stress tensor. This non-perturbative approxima- 
tion extrapolates the Lagrangian linear model by the re- 
placement Ax — > Ax+p- 
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There are other models that are not covered by this 
list, notably higher-order perturbative approximations. 
For example, in the Lagrangian perturbation framework 
that contains the Zel'dovich approximation as a special 
first-order solution in a subclass (Buchert 1989, 1992), so- 
lutions have been derived for dust up to the fourth order. 
Special classes of second-order solutions and their inves- 
tigation may be found in Bouchet et al. (1992, 1995) and 
Buchert & Ehlers (1993). To the third order Moutarde 
et al. (1991) gave a special solution, Buchert (1994) and 
Catelan (1995) investigated a large class. Vanselow (1996) 
derived second- and third-order solutions for some more 
general cases, as well as a class of fourth-order solutions. 
By including isotropic stresses, Lagrangian perturbation 
solutions were given to first order (Adler & Buchert 1999), 
to second order (Morita & Tatekawa 2001; Tatekawa et al. 
2002), and to third order (Tatekawa 2005a,b). All of these 
higher-order models can be regarded as approximating the 
residual vector field TZ. However, in all of these models, 
the peculiar-gravitational field strength remains smooth 
while crossing high-density regions, which (as becomes 
clear below) points to a shortcoming of the perturbative 
calculations. The difficulty in deriving a closed differen- 
tial equation for the displacement field P in the highly 
non-linear regime lies precisely in the term 7Z, which is a 
non-local (in space and time) and a non-linear functional 
of P; approximating 7Z perturbatively hides an important 
effect, as we shall see in detail below. 

Finding an extension into the non-perturbative (both 
Eulcrian and Lagrangian) regimes is an involved math- 
ematical task. Notwithstanding, it should be attempted. 
For example, even the Lagrangian perturbation approach 
falls short capturing the action of multi-stream forces. 
This was demonstrated by comparing the statistical prop- 
erties of second- and third-order Lagrangian perturbation 
solutions with results of numerical simulations (Tatekawa 
2004a, b). This shortcoming calls for a non-perturbative 
generalization. 

A proper understanding of the extrapolation of the 
linear relationship (|10f> between w and P into the non- 
perturbative regime requires that the residual term 7Z. in 
Eq. (|6a[) be analyzed in general. To this end we exploit 
the fact that Eq. l|6a[) holds independently of whether we 
are talking about dust or a general dispersion-supported 
system. It is possible to find exact integrals of the general 
equation l|6a[) . which will also help us to understand the 
quality of the relationship (fTTTf) - lying at the basis of most 
currently known models - without involving higher-order 
perturbation analysis. 

Finally, although virialized states can be understood 
through the tensor virial theorem, eventually including 
surface terms to account for the non-isolated state of 
gravitational systems, there are signatures possibly im- 
printed onto the phase space distribution during the non- 
perturbative regime. This distribution may have a "re- 
laxed" global shape, but its internal structure will prob- 
ably appear structured, rather than completely smooth, 



e.g. as a result of a hierarchy of embedded (smoothed) 
caustics (Ed Bertschinger, priv. comm.). 

3. Exact integral for the peculiar-gravitational 
field strength 

In this section we look at the general equations and com- 
bine them in order to obtain a transport equation for the 
peculiar-gravitational field strength. We then integrate 
the transport equation exactly with the help of a (suffi- 
cient) restricting condition. With this result we are able to 
calculate the gravitational field strength for a given family 
of trajectories, and to solve the longitudinal field equation 
V q • w = —A-kGqhclS in general. We later demonstrate 
that this integral can be exploited to obtain a powerful 
approximation for the non-perturbative regime of struc- 
ture formation. 

3.1. Transport equation for the peculiar-gravitational 
field strength 

We start with Eq. @ and employ the vector identity 

(u ■ V q )w - u(V q ■ w) = 
[(w • V q )u- w(V q • u)] + V q x (w x u) (13) 

to find the modified evolution equation: 



1 

O L 



— w + 2Hw — AttGqhvl = 
dt 

( W • V q ) U - W ( V q • U ) + V q X f 



(14) 



with a new vector potential T := T + w x u. Computing 
@2t ( w /fiO ' an d using the continuity equation (|laf) . we 
arrive at the following equation that we may call transport 
equation for w: 



w 

— • V, 



d 








(t)-*< 




dt 




V Q J 


u 


A ^ U 


_JLi 


- -+ 


- 4ttGq h - H 




a 


Q 


ag 



(15) 



This is the key-equation of the present work. 



3.2. Integrating the transport equation 

For the case V q x T = 0, we can find an exact integral to 
the above transport equation along comoving trajectories 
q = F(X, t) as follows. 

Motivated by a recent investigation of an exact 
Lagrangian integral for the gravitational field strength g in 
Newtonian gravity (Buchert 2006a), we make the follow- 
ing ansatz (which generalizes the ansatz for the case of a 
non-vanishing background source in the above-mentioned 
work) : 



w - CF 



a(k-V )F ; C=C(*) 



k = k(X) , 
(16) 
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with the nabla-operator with respect to Lagrangian coor- 
dinates Vo- Performing the total time-derivative of this 
ansatz along the integral curves F = X + P of the scaled 
peculiar-velocity field, i.e., = ^u(F, t), using a fur- 
ther identity (e.g., Serrin 1959), 

(k-V ) jF = [(k-V )F].V q " , (17) 

and once more applying the identity (|13fl to the fields u 
and F, 



V q - [F(V q -u)-(F-V q )u] = 
[u(V q ] 



(u-V q )F] = V q -(d-l)u 



(18) 



(with d denoting the dimension of space), we obtain an 
equation that we compare with the transport equation. 
We are left with the following conditions on the unknown 
function £: 



F 

Q 



(( + 2HC 







C 



AttGqho 







For non-vanishing F and u, the above two equations for 
£ are equivalent by virtue of qh = SHqh, and we have 
determined the unkown function: 

AttGqh a 



c 



(19) 



Now, we can write down an exact integral w 1 for the 
peculiar-gravitational field-strength. We also replace the 
density by its exact Lagrangian integral, 



e (x,t) 



gPMo) 
J 



Qh- 



l + 5(X,t ) 
Jf 



QH 



QH(to) 



(20) 

where J := det(/ji/.) denotes the Jacobian of the transfor- 
mation from x = f(X,t) to X, and Jp := det(F^ k ) = 
Ja~ 3 denotes the Jacobian of the transformation from 
q = F(X, t) to X with the comoving Lagrangian defor- 
mation gradient (F^) 5 ). The result reads 



K = 



w J = 



= W 



( K ■ V ) F ^Ggna ] 
o?Jf d 

47rGgjj(t ) - 



(21) 



— - — — ,, , -X is the integration constant, 

g(t ) d o ' 

which was determined by evaluating the integral at t = to . 
The integral 121J1 explicitly depends on the dimension d of 
the continuum. 

The above integral is quasi-local, i.e. it locally repre- 
sents the peculiar-gravitational field strength through a 
functional of the comoving trajectory field F, while initial 
data are constructed non-locally according to the struc- 
ture of the theory. 

A strong property of the integral l|21|) is that it solves 
the longitudinal field equation V q • w = —AttGqhoS in 
general. (We can also find such an integral for the remain- 
ing field equations V q x w = 0, see Sect. 5; a general 
integral is expected to include non-local terms.) 

5 We denote a spatial derivative with respect to Lagrangian 
coordinates by a vertical slash | that commutes with the 
Lagrangian time-derivative. 



3.3. Lagrangian framework, proof of the exact integral 
and its transformation properties 

Let us recall the basic equations of a Lagrangian descrip- 
tion of fluid motion. As shown in Buchert (2006a), the pre- 
sented integral (|21|l can also be obtained by simply trans- 
forming the result obtained in Newtonian gravity. This is 
expected for exact equations and expressions. The explicit 
derivation in this paper shows that the approximation of 
a vanishing curl of the vector potential T apparently does 
not impair this property; in general, approximations to 
the system of equations in Newtonian gravity do not carry 
over to corresponding approximations in Newtonian cos- 
mology by a simple transformation. 

In Newtonian gravity the Lagrangian description is 
based on integral curves x = f(X, t) of the full velocity 
field v(x, t): 



§ = «'•<> 



f(X,to) =:X 



(22) 



By introducing this family of trajectories, we can express 
all Eulerian fields (e.g., the velocity v, the acceleration g, 
the density q, and the vorticity u> := ^V x x v) in terms 
of the field of trajectories x = f (X, t) as follows: 



v = f(X,t) ; g = f(X,t) 



Q = 



00 
J 



■ Vpf 

J 



(23) 



(24) 



with the Jacobian of the transformation from Eulerian to 
Lagrangian coordinates J := det(/jife(X, t)) > 0, qq := 
e(X,io), w := w(X,i ). 

Equation l|24(l lists the known Lagrangian integrals for 
the density q and the vorticity us of the Euler-Newton 
system; i.e., they represent a Eulerian field as a functional 
of f . To transform those fields back to Eulerian space, we 
need the transformation f to be invcrtiblc: i.e. J > 0, 
defining regular solutions (for more details the reader may 
consult the review by Ehlers & Buchert 1997). 

The Eulerian field equations are transformed into a 
system of Lagrangian equations by virtue of the following 
transformation of the field strength gradient: 

~Q^9i = 9i\k^-hk= 7pf e Hm^jpq9i\kfp\tfq\m ■ (25) 



2J 



The gradient of the inverse transformation from 
Lagrangian to Eulerian coordinates, h = f , was ex- 
pressed in terms of f through the algebraic relationship 



dx i 



'^ij — dd{J%j)J — j€ik££jmnfm\kfn\i 



1 

2J" 



(26) 

For the field equations we obtain with (|25() the following 
set of four Lagrangian equations (Buchert & Gotz 1987 
for A = 0, and Buchert 1989 for A ^ 0) (i,j,k = 1,2,3 
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with cyclic ordering; summation over repeated indices is 
understood) : 



1 9(g a J b J c ) 
2 £abc 8(X 1 ,X 2 ,X 3 ) 



e pq[j 



d(g t ]J P Jg) 
d(x 1 ,x 2 ,x 3 ] 



-AJ = -4vrG£»o(X) 







(27) 
(28) 



In the case of a dust continuum, the above equations form 
a set of four evolution equations, the Lagrange-Newton 
system of equations, by virtue of g = f . Alternative forms 
of these equations may be found in Buchert (1996) and 
Ehlers k Buchert (1997). 

The above equations not only hold for a continuum of 
dust; they simply represent a transformation of the field 
equations, which also hold in the more general setting 
discussed in Sect. 2. The difference comes in when we 
represent the field strength g in terms of the trajectory 
field; e.g. including isotropic pressure forces the resulting 
Lagrangian evolution equations are investigated in Adler 
& Buchert (1999). 

Proposition 2 in Buchert (2006a) gives the following 
integral of the first Lagrangian equation (including its 
proof): 



A 



g /A = (C • Vo)f J 1 , C:=G 



(29) 



Since this result is exact, we are entitled to apply the 
transformation to comoving coordinates and peculiar- 
fields to this result. Since we can adopt the same 
Lagrangian coordinates X in both cases in view of x = 
a(t)q, a(to) = 1 and since the background field strength 
evolves as gn = aq, we make the ansatz g = aF+w(X, t), 
G = a(to)X + W(X) and use Friedmann's equation 
3d = A — A-kGoqh to obtain the integral w 7 in Eq. (|21|l 6 . 

4. The exact integral and its implications 

In order to learn more about the implications of the in- 
tegral Q]} , 

we now recover known model equations in 
cosmology (e.g. the Lagrangian perturbation scheme; this 
should be possible for the longitudinal parts, since the in- 
tegral (|21|l is general in this case). Then, we elaborate on 
possible applications. 

4.1. Recovering known cosmological models 

We first write the general expression for w in terms of P. 
In the case of a dust continuum, w = a(P -I- 2HP) (from 
Euler's equation for dust w = u + Hu with u = aP), we 
rewrite the integral l|21|) in terms of P (we drop the index 
/ to ease the notation, consider the case d = 3 and divide 
by a): 

P + 2HP = (K-Vo)(X + P) + 4 ! rG g H 
o\l F 3 



6 For an independent proof we could instead insert 121H into 
the transformed Lagrangian equations given in Appendix A of 
Buchert (1989). 



with the following general expresssion for the Jacobian 
J F = a~ 3 J: 



J F = l + I(P m ) + II(P i{k ) + III(P tlk ) 



(31) 



where /, II and 777 denote the principal scalar invariants 
of the tensor in brackets. 

Let us first expand the expression l/J F with the 
Jacobian J23} to first order, 1/J F = 1/(1 + V • P + ...) ~ 



1 - V • P. We insert K =: W 



47rGg g (t ) 



X, drop terms 



that are quadratic in the fields P and W, and replace the 
term (X • Vo)P — X(Vo • P) by the (up to a transversal 
part equivalent) term (P-V )X-P(V -X) = -2P. This 
approximated equation then reads: 



pW + 2HP (1) - 4irGg H P W = 



W 



(32) 



Note that the factor 3 disappeared from the term in front 
of P' 1 '. This is the equation for longitudinal first-order 
Lagrangian perturbations (Buchert 1989, 1992). For ini- 
tial quasi -homogeneity, which is sometimes assumed, cf. 
Footnote 4, the right-hand side of this equation drops. 

Expanding l/J F to higher orders soon yields messy 
expressions; e.g. by keeping only terms quadratic in PW 
and W, we deal with an approximate equation of the form 



p(2) + 2H p{2) _ A7rGgH p(2) = W 

a 6 

+ 4 f W • V P (1) - WV • P (1) 
a 6 V 



(33) 



The last term in the above equation illustrates that we 
recover source terms for the second-order perturbation 
solution P( 2 ', which correspond to the local parts (as in- 
troduced in Buchert 1993, Buchert 1994); the full terms in- 
cluding the non-local parts are solutions of Poisson equa- 
tions with sources given by the divergence of expressions 
of this type. 

By expressing w through more general Euler equa- 
tions that include presssure terms, velocity dispersion, 
or deviations from mean field gravity (cf. Buchert & 
Dominguez 2005), we can also expand the integral in or- 
der to recover known approximations, e.g., the first-order 
Lagrangian equation for longitudinal perturbations in a 
medium supported by isotropic pressure forces (Adler & 
Buchert 1999). 

4.2. Iteration approach 

A powerful possibility of applying the exact integral fo- 
cusses on an iterative definition of the trajectories or of 
the displacement field. This iterative view exploits the fact 
that we know the field strength exactly (and in general 
with respect to the longitudinal field equation) along any 
trajectory field that we could imagine. 

We now exemplify the iteration procedure for a con- 
tinuum of dust. Note that the integral (|21l) can be viewed 
as a set of three partial differential equations for the tra- 
jectory field after inserting the general expression for the 
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peculiar-acceleration field, e.g. for dust: w = a(F + 2HF). 
Holding the Lagrangian coordinates fixed, i.e. for a single 
trajectory, we are dealing with ordinary differential equa- 
tions. 

The way we define the iteration scheme enjoys some 
freedom (non-unique definition). We argue below for the 
following choice of writing the integral pifl as an iteration 
scheme: 



a 



(K-Vq)fW 8ttGqh j 
a 3j|? ] 3 



(34) 

with K given in and iw'^ 1 ' = F^ 1 ! + 2H¥^ n+1 \ 
denoting the iteration steps in brackets to avoid confusion 
with the perturbation index used earlier. In this way we 
have split the linear term 4 "° te F for d = 3 into AnGg H - 
^AttGqh, the first term we consider as belonging to [n+ 1] 
and the second to [n] . The reason for this choice becomes 
obvious after the remarks below. 

For the displacement field P = F — X, we obtain the 
following iteration scheme: 



p[n+i] + 2iJp[«+ 1 l - 4nGg H P ln+1 



A-kGqh 



X - 2P l; 



(K-V )(X + PN) 



i 3 det 



(Si, 



i\k 



(35) 



Iteration consists in the strategy of feeding in a comoving 
displacement field P^ on the right-hand side of Eq. (|35J) 
in order to obtain another displacement field p[™ +1 l by 
solving second-order, ordinary differential equations for 
each fluid element X. The above choice of the iteration 
scheme is motivated by three arguments. First, starting 
the iteration l|34l) with the trivial trajectory field F^ = X, 
which corresponds to a straight Hubble expansion, we ob- 
tain the equation governing the longitudinal first-order 
Lagrangian approximation exactly, F^ = X + pW, with 
pW obeying Eq. Ij32(l for P' 1 ). Second, the choice is sup- 
ported by the fact that F^ (as a special exact solution 
of the Lagrange-Newton system {127128$ -) also produces 
a special exact solution FN of the same system (Buchert 
1989). Third, to all orders in Lagrangian perturbation the- 
ory, the longitudinal differential operators for p( n+1 ) have 
all the form of the differential operator for pW, Eq. I|32l) . 
with sources involving the lower-order perturbation solu- 
tions (see Ehlers & Buchert 1997, Sect. 3.2, for the gen- 
eral perturbation and solution schemes); compare with the 
left-hand side of Eq. (J2SJ. 

However, this does not mean that further iteration 
eventually produces further exact solutions, and also that 
another choice could not perform as well as the above 
choice. In order to illustrate the iteration scheme further, 
we compute the equation for the second iterate P^ 2 ' . While 
the first iterate (insert P^ =0 into Eq. 1)3501 obeys 



(36) 



P [i] + 2HPW - AttGqhPW = W f 



the second iterate may be found by inserting only a sub- 
class of the general solution to Eq. (|3l)|) that corresponds to 



Zel'dovich's approximation (Zel'dovich 1970, 1973). After 
factorizing pM into time-dependent functions and vector- 
functions of initial data, the time-dependent solutions 
consist of two homogeneous and one particular solution 
of the following equation (Buchert 1992): 

£(t) + 2J7(t)£(t)-47rG0 ff (t)(£(t)+l) = , (37a) 

where H(t) = ~ and gn(t) = QH{to)a~ 3 have to be ex- 
pressed through solutions of Friedmann's differential equa- 
tion: 



inGgn 



(37b) 

a o a z 3 

Explicit forms of the functions £(i) including a cosmo- 
logical constant can be found in (Bildhauer et al. 1992; 
see also the supplement by Chernin et al. 2003). We 
choose the restricted initial data set for which U oc W 
(U(X) := u(q,io)) (for more details see: Buchert 1992). 
Inserting this subclass of first-order solutions (Zel'dovich's 
approximation), Pl lz I := 6(t)Vo*(X), with the growing 
mode solution =: b(t) of Eq. (|37a|l . we can find the 
second iterate pl 2Z l corresponding to this restricted choice 
as follows (we also use expression H3H ): 

p[2Z] + 2HV [2Z] - 4irGg H P l2Z] = 
AitGqh 



■(X-26(t)V tf(X)) 4 

o 

(K- V )(X + &(t)V fr(X)) 

[ i + &/(*| ifc ) + 6 2 //(*| 4fc ) + bnii{*\ ik ) ] 



(38) 



This equation defines a set of ordinary differential equa- 
tions parametrized by X. Its analytical solution (the sub- 
ject of a forthcoming work) provides a test case for a nu- 
merical iteration scheme. For this purpose we now give 
some useful technicalities. 



4.3. Numerical implementation 

For the purpose of numerically implementing the iteration 
scheme l|35|l . we may rescale the dependent and indepen- 
dent variables as follows. Following Shandarin (1980) ex- 
cept that we refer all quantities with subscript "0" to the 
initial time we introduce the dimensionless and appro- 
priately scaled variables (Buchert 1989, App.A) 

q:=q/<7o I u := (u/u )a ; w := (w/w )a 3 , 

u = qo/t ; w = qo/tl ; g := (g/ QH(to))a 3 , (39a) 

and a conformal transformation of the time-variable: 

dT:=±^- , (39b) 
to a 2 (t) 

which is negative and tends to — oo at the Big-Bang singu- 
larity. With the help of (|39b|l . we can write the solutions 
of Friedmann's differential equation i|37b(l and the mass 
density parameter for the cases A = in simple forms: 



a(T) 



T 2 + k 



fi r ' 



3H 2 



rp2 



(39c) 



7 For the cases A / we have to employ other strategies, 
e.g.: Bildhauer et al. 1992. 
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with k = 0,±1. For an Einstein-de Sitter model (k = 
0), we have To = —3, and for the other models To = 

w-rj- Below we also use the relations 4ttGqh (to) — 



^H 2 and A7rGg H (to)t 2 = 6/(T 2 + k)\ 
(for further details see: Buchert 1989, App.A). 
For the scaled displacement field, 



P:=P 



X := — ; V = <jo V 
<7o 



(39d) 



Eq. 1(35(1 can be rewritten as a set of first-order (for each 
X ordinary) differential equations: 



dT 



T 2 + k d _, 



_,-,["+!] _ Q p[n+l] 

2 dT U ^ 



X - 2P [ " 



K ■ V X + P 



det S. 



(40) 



with K := i(T 2 + fc)W - X. The differential operator in 
the second equation simplifies further by integrating along 
the time r with dr := 2a(T)/(T 2 + k) dT. Initial data for 
W follow from solving the initial Poisson equation for the 
density contrast 8 = g — 1: 



fc)W = -3 <S(X) ; V x W = . (41) 



Initial data for U could be specified by the special choice 
U = W. 

As a first test of a numerical scheme one chooses 
PM = 0, so that PW should be identical to Zel'dovich's 
approximation (Zel'dovich 1970, 1973). A second test will 
be provided by the analytical solution for the second iter- 
ate, as mentioned above. A third test should measure the 
(artificial) vorticity of w as introduced by the integral for 
generic initial data; here, one can compare with the exact 
expression (|44|l below for the transverse part of the field 
strength. 

All peculiar-fields appearing in the calculations must 
be periodic on the largest scale to assure that the aver- 
age model is Friedmannian. This construction is necessary 
for the uniqueness of a Newtonian solution and implies a 
globally vanishing 'backreaction' (see: Buchert & Ehlers 
1997). 



peculiar-velocity, which is a very good assumption in the 
weakly non-linear regime, but certainly fails in a highly 
non-linear situation. Consequences of the above remarks 
have been investigated in detail in (Buchert & Dommguez 
2005). 

We have argued that relationship (fTUf> needs gener- 
alization for the understanding of the non-perturbative 
regime of cosmic structure formation. Integrating the 
transport equation for the peculiar-gravitational field 
strength, Eq. 1(15(1 , we obtained an exact integral, general- 
izing relationship l|10fl ■ It was demonstrated that this inte- 
gral can be employed to define an iteration scheme that al- 
lows us to obtain the peculiar-gravitational field strength 
for any given family of trajectories, which solves the 
Lagrangian evolution equation corresponding to V q • w = 
—AitGqhciS in general. This property, together with the 
experience of the good performance of Lagrangian pertur- 
bation schemes based on the longitudinal part only (e.g., 
Buchert et al. 1997), supports the expectation that the 
integral and its corresponding iteration scheme provide a 
powerful approximation for the non-perturbative regime 
of structure formation. 

The qualitative difference from a Lagrangian pertur- 
bation analysis is due to the fact that the former lacks 
the important leading term proportional to the density in 
1)21(1. Lagrangian perturbation solutions predict a smooth 
gravitational field strength when crossing caustics in the 
density field for all orders in perturbation theory (com- 
pare the general perturbation and solution schemes given 
by Ehlers & Buchert 1997). Therefore, we are entitled to 
consider the integral (|21|l as a genuinly non-perturbative 
result. It shows that a blow-up of the field strength at 
caustics is a generic property of the gravitational collapse. 
Counterarguments based on exact solutions with sym- 
metry do not apply to the generic situation. Consider, 
as an example, plane-symmetric motions on a three- 
dimensional homogeneous-isotropic background. In that 
case the general exact solution (at the same time a solu- 
tion of Eq. ((32(1 ) does not predict a singular field strength: 
specifying the integral ((21(1 to plane-symmetric motion, 
we infer from the expression proportional to the density, 



((W 



A-irGg H {to) - 



X) • V (X + P) 



q 2 Jf 



(42) 



5. Summary and prospects 

Current analytical models (excluding higher-order pertur- 
bative corrections) basically follow from the simple as- 
sumption that the peculiar-gravitational field strength 
is proportional to the displacement field, Equation l(10|) : 
w = + 4:TtGqhciP- Eq. ((12(1 furnishes the up to date 
most general model equation based on this assumption. 
More frequently employed models, like the celebrated 
Zel'dovich approximation and the adhesion approxima- 
tion, even imply the tighter restriction of proportional- 
ity for the peculiar-gravitational field strength to the 



with J F = l + P xll , (43) 

that the Jacobian cancels the directional derivative term 
exactly, and w; 1 plane = AttGqhciPi + W\/a 2 does not blow 
up at the caustic where 1+ Pi|i =0. 

One of the most promising application fields of the 
integral and its iteration scheme could be the following. 
Since 1(211) provides w as a local functional of the dis- 
placement field P, it may substantially enhance the power 
of reconstruction methods (Croft & Gaztanaga 1997; 
Susperregi & Buchert 1997; Courteau & Willick 2000, 
Courteau & Dekel 2001; Brenier et al. 2003, Mohayaee 
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et al. 2003). As mentioned above, the integral (|2*T)l gener- 
alizes the assumption l|10[) . which itself lies at the basis of 
most analytical models (as summarized in this paper) and 
which already furnishes a more general ansatz for recon- 
struction methods, compared to those that are currently 
implemented, e.g. those that assume the proportionality 
of w and u. 

To propose a non-perturbative approximation based 
on the integral (|21|1 is also motivated by simplicity, be- 
sides the advantage of a local expression. On theoreti- 
cal grounds, however, there are limitations to the inte- 
gral J2U, because it will produce artificial vorticity of 
w, since a generic approximation based on (|21l) will not 
satisfy V q x w = 0. Numerical work has to confirm the 
expectation that the longitudinal part of the field equa- 
tions provides the dominant contribution. The transverse 
contribution and the contribution from non-local terms 
can also be quantified, based on results of a forthcom- 
ing work on transported differential forms associated with 
the gravitational field strength; we here already give the 
exact Lagrangian integral that solves the transverse field 
equation in general (Buchert 2006b): 



r trans verse 



1 

aJi 



sub 



(Wx [(V xF)-Vo P])*), 
(44) 

where J™ b denotes the subdeterminants of (i^ifc)- 

Since a general integral, if it exists, would include non- 
local parts, and thus would require the solution of elliptic 
boundary value problems at all times, it may be of limited 
practical use. 

Acknowledgements. This work was supported by the 
"Sonderforscmmgsbereich SFB 375 fur Astro-Teilchenphysik 
der Deutschen Forschungsgemeinschaft" . Thanks go to 
Stephane Colombi for stimulating discussions and thoughts 
on the numerical implementation of the iteration scheme. 

References 

Adler, S., & Buchert, T. 1999, A&A, 343, 317 

Bagla, J.S., & Padmanabhan, T. 1997, Pramana, 49, 161— 

192 | |astro-ph/0411730| | 
Bcrnardeau, F., Colombi, S., Gaztahaga, E., & 

Scoccimarro, R. 2002, Phys. Rep., 367, 1 
Bertschinger, E. 1998, ARA&A, 36, 599 
Bildhauer, S., Buchert, T., & Kasai, M. 1992, A&A, 263, 

23 

Binncy, J., & Tremaine, S. 1987, Galactic Dynamics, 

Princeton University Press 
Bouchet, F.R., Juszkiewicz, R., Colombi, S. and Pellat, R. 

1992, ApJ, 394, L5 
Bouchet, F.R., Colombi, S., Hivon, E. and Juszkiewicz, R. 

1995, A&A, 296, 575 
Brcnier, Y., Frisch, U., Henon, M., Loeper, G., Matarrese, 

S., Mohayaee, R., & Sobolevskii, A. 2003, MNRAS, 

346, 501 
Buchert, T. 1989, A&A, 223, 9 
Buchert, T. 1992, MNRAS, 254, 729 
Buchert, T. 1993, A&A, 267, L51 



Buchert, T. 1994, MNRAS, 267, 811 

Buchert, T. 1996 In Proc. IOP 'Enrico Fermi', Course 

CXXXII (Dark Matter in the Universe), Varenna 1995, 

eds. S. Bonometto, J. Primack, & A. Provenzale, IOS 

Press Amsterdam, pp. 543-564 
Buchert, T. 2006a, Phys. Lett. A, 354, 8 
Buchert, T. 2006b, work in preparation 
Buchert, T., & Dommguez, A. 1998, A&A, 335, 395 
Buchert, T., & Dommguez, A. 2005, A&A, 438, 443 
Buchert, T., & Ehlers, J. 1993, MNRAS, 264, 375 
Buchert, T., & Ehlers, J. 1997, A&A, 320, 1 
Buchert, T., & Gotz, G. 1987, J. Math. Phys., 28, 2714 
Buchert, T., Karakatsanis, G., Klaffl, R., & Schiller, P. 

1997, A&A, 318, 1 
Catelan, P. 1995, MNRAS, 276, 115 

Chernin, A.D., Nagirner, D.I., & Starikova, S.V. 2003, 
A&A, 399, 19 

Courteau, S. & Dekel, A. 2001, in: Astrophysical Ages and 

Times Scales, ASP Conference Series Vol. 245., T. v. 

Hippel, C. Simpson, N. Manset (eds.) 
Courteau, S. & Willick, J. 2000 (eds.), Cosmic Flows 

Workshop, ASP Conference Series, Vol. 201. 
Croft, R.A.C., & Gaztanaga, E. 1997, MNRAS, 285, 793 
Ehlers, J., & Buchert, T. 1997, Gen. Rel. Grav., 29, 733 
Gabrielli, A., Baertschiger, T., Joyce, M., Marcos, B. & 

Sylos-Labini, F. 2006; |cond-m at/060312^ 
Gurbatov, S.N., Saichev, A.I., & Shandarin, S.F. 1989, 

MNRAS, 236, 385 
Ma, C.-P, & Bertschinger, E. 2004, Ap.J. 612, 28 
Morita, M., & Tatekawa, T. 2001, MNRAS, 328, 815 
Mohayaee, R., Frisch, U., Matarrese, S., & Sobolevskii, A 

2003, A&A, 406, 393 
Moutarde, F., Alimi, J.-M., Bouchet, F.R., Pellat, R., & 

Ramani, A. 1991, Ap.J. 382, 377 
Peebles, P.J.E. 1980, The Large-scale Structure of the 

Universe, Princeton Univ. Press 
Sahni, V., & Coles, P. 1995, Phys. Rep. , 262, 1 
Shandarin S.F. 1980, Astrofizika 16, 769; 1981: 

Astrophysics 16, 439 
Serrin, I. 1959 In Encyclopedia of Physics, Vol. III. 1 , ed. 

S. Flugge, Springer Berlin 
Shirokov, A., & Bertschinger, E. 2006, submitted to 

Ap. J.Suppl.; |astro-ph/0505087| 
Stcinmctz, M., 2003, in: Hubble's Science Legacy: 

Future Optical/ Ultraviolet Astronomy from Space, 

ASP Conference Proceedings, K.R. Sembach et al. 

(eds.) Vol. 291, p.237 
Susperregi, M., & Buchert, T. 1997, A&A 323, 295 
Tatekawa, T., Suda, M., Maeda, K., Morita, M., & Anzai, 

H. 2002, Phys. Rev. D66, 064014 
Tatekawa, T. 2004a, Phys. Rev. D69, 084020 
Tatekawa, T. 2004b, Phys. Rev. D70, 064010 
Tatekawa, T. 2005a, Phys. Rev. D71, 044024 
Tatekawa, T. 2005b, Phys. Rev. D72, 024005 
Vanselow, M. 1995, Diploma Thesis, Ludwig-Maximi- 

lians-Universitat Mimchcn (in German). 
Zel'dovich, Ya.B. 1970, A&A, 5, 84 
Zel'dovich, Ya.B. 1973, Astrophysics, 6, 164 



